{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eddy detection and filter\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\n\nfrom matplotlib import pyplot as plt\nfrom numpy import arange\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def start_axes(title):\n    fig = plt.figure(figsize=(13, 5))\n    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n    ax.set_aspect(\"equal\")\n    ax.set_title(title, weight=\"bold\")\n    return ax\n\n\ndef update_axes(ax, mappable=None):\n    ax.grid()\n    if mappable:\n        plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load Input grid, ADT is used to detect eddies.\nAdd a new filed to store the high-pass filtered ADT\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = RegularGridDataset(\n    data.get_demo_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"),\n    \"longitude\",\n    \"latitude\",\n)\ng.add_uv(\"adt\")\ng.copy(\"adt\", \"adt_high\")\nwavelength = 800\ng.bessel_high_filter(\"adt_high\", wavelength)\ndate = datetime(2016, 5, 15)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the detection for the total grid and the filtered grid\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a_filtered, c_filtered = g.eddy_identification(\"adt_high\", \"u\", \"v\", date, 0.002)\nmerge_f = a_filtered.merge(c_filtered)\na_tot, c_tot = g.eddy_identification(\"adt\", \"u\", \"v\", date, 0.002)\nmerge_t = a_tot.merge(c_tot)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display the two detections\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Eddies detected over ADT\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nmerge_f.display(\n    ax,\n    lw=0.75,\n    label=\"Eddies in the filtered grid ({nb_obs} eddies)\",\n    ref=-10,\n    color=\"k\",\n)\nmerge_t.display(\n    ax, lw=0.75, label=\"Eddies without filter ({nb_obs} eddies)\", ref=-10, color=\"r\"\n)\nax.legend()\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Amplitude and Speed Radius distributions\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 5))\nax_a = fig.add_subplot(121, xlabel=\"Amplitude (cm)\")\nax_r = fig.add_subplot(122, xlabel=\"Speed Radius (km)\")\nax_a.hist(\n    merge_f.amplitude * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddies in the filtered grid\",\n    histtype=\"step\",\n)\nax_a.hist(\n    merge_t.amplitude * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddies without filter\",\n    histtype=\"step\",\n)\nax_a.set_xlim(0, 10)\nax_r.hist(merge_f.radius_s / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.hist(merge_t.radius_s / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.set_xlim(0, 100)\nax_a.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Match detection and compare\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_, j_, c = merge_f.match(merge_t, cmin=0.1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Where are the lonely eddies?\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kwargs_f = dict(lw=1.5, label=\"Lonely eddies in the filtered grid\", ref=-10, color=\"k\")\nkwargs_t = dict(lw=1.5, label=\"Lonely eddies without filter\", ref=-10, color=\"r\")\nax = start_axes(\"Eddies with no match, over filtered ADT\")\nmappable = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nmerge_f.index(i_, reverse=True).display(ax, **kwargs_f)\nmerge_t.index(j_, reverse=True).display(ax, **kwargs_t)\nax.legend()\nupdate_axes(ax, mappable)\n\nax = start_axes(\"Eddies with no match, over filtered ADT (zoom)\")\nax.set_xlim(25, 36), ax.set_ylim(31, 35.25)\nmappable = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nu, v = g.grid(\"u\").T, g.grid(\"v\").T\nax.quiver(g.x_c, g.y_c, u, v, scale=10, pivot=\"mid\", color=\"gray\")\nmerge_f.index(i_, reverse=True).display(ax, **kwargs_f)\nmerge_t.index(j_, reverse=True).display(ax, **kwargs_t)\nax.legend()\nupdate_axes(ax, mappable)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot ({i_.shape[0]} matches)\", weight=\"bold\")\n\nfor i, (label, field, factor, stop) in enumerate(\n    (\n        (\"Speed radius (km)\", \"radius_s\", 0.001, 80),\n        (\"Effective radius (km)\", \"radius_e\", 0.001, 120),\n        (\"Amplitude (cm)\", \"amplitude\", 100, 25),\n        (\"Maximum Speed (cm/s)\", \"speed_average\", 100, 25),\n    )\n):\n    ax = fig.add_subplot(\n        2, 2, i + 1, xlabel=\"Filtered grid\", ylabel=\"Without filter\", title=label\n    )\n    ax.plot(merge_f[field][i_] * factor, merge_t[field][j_] * factor, \".\")\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.plot((0, 1000), (0, 1000), \"r\")\n    ax.set_xlim(0, stop), ax.set_ylim(0, stop)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}